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ABSTRACT 

We present a method to constrain galaxy parameters directly from three-dimensional data cubes. 

The algorithm compares directly the data with a parametric model mapped in x , y, A coordinates. It 
uses the spectral line-spread function (LSF) and the spatial point-spread function (PSF) to generate 
a 3-dimensional kernel whose characteristics are instrument specific or user generated. The algorithm 
returns the intrinsic modeled properties along with both an ‘intrinsic’ model data cube and the mod¬ 
eled galaxy convolved with the 3D-kernel. The algorithm uses a Markov Chain Monte Carlo (MCMC) 
approach with a nontraditional proposal distribution in order to efficiently probe the parameter space. 

We demonstrate the robustness of the algorithm using 1728 mock galaxies and galaxies generated from 
hydrodynamical simulations in various seeing conditions from 0'/6 to 1/2. We find that the algorithm 
can recover the morphological parameters (inclination, position angle) to within 10% and the kine¬ 
matic parameters (maximum rotation velocity) to within 20%, irrespectively of the PSF in seeing 
(up to 1/2) provided that the maximum signal-to-noise ratio (SNR) is greater than ~ 3 pixel -1 and 
that ratio of galaxy half-light radius to seeing radius is greater than about 1.5. One can use such an 
algorithm to constrain simultaneously the kinematics and morphological parameters of (nonmerging) 
galaxies observed in nonoptimal seeing conditions. The algorithm can also be used on adaptive-optics 
(AO) data or on high-quality, high-S/N data to look for nonaxisymmetric structures in the residuals. 

Subject headings: methods: data analysis — methods: numerical — techniques: imaging spectroscopy 


1. INTRODUCTION 

Thanks to several studies using optical or near-infrared 
(NIR) integral field unit (IFU) spectroscopy of Ha emis¬ 
sion from local and high-redshift (z > 1) galaxies 


(Forster Schreiber et al. 2006] Law et al. 

2007 van 

Starkenburg et al. 20U8{ Cresci et al. 
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rnation has changed significantly in the past decade. For 
instance, these surveys have shown that a significant sub¬ 
set of high-redshift galaxies have a disklike morphology 
and show organized rotation, with regular velocity field s. 


In contrast to low-redshift studies (e.g. Bacon et al. 
2001 Cappellari et al. 20111, high-redshift (1 z IS ^) 


galaxies are observed at a spatial resolution that is sever- 
aly limited by the seeing conditions owing to their small 
apparent angular sizes. In order to overcome the low spa- 


are often required (Law et al. 

2007 

2008 

2011 Wright et al. 2007 

2009 


2009 Genzel et al. 


2009). However, observa- 


rnental costs, and add strong observational constraints 
such as the additional exposure times required to com- 
pensate for the loss in surface brightness (SB) sensitivity 
(Law et al7||2006 1. Indeed, the SB limit for AO observa- 


tions taken on smaller pixels is higher, leaving the cur¬ 
rent state-of-the art observations to the objects with the 
highest SBs. 
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Given these challenges and the advancements in mul¬ 
tiplexing IFU observations with the Very Large Tele- 
scope (VLT) second-g eneration instruments like KMOS 
(Sharpies et al.| 2006 ) and the Multi-unit Spectrograph 
Explorer (MUSE; Bacon et al. 2006 2015), it is impor¬ 
tant to have tools that can give robust estimates on 
the galaxy physical properties. In particular, KMOS 
will bring large statistically significant samples of high- 
redshift galaxies as it can observe 24 galaxies at a time, 
but this facility will always lack an AO unit. This could 
potentially be a serious limitation since the robustness 
of the derived kinematic parameters may depend on the 
quality of the atmospheric conditions (seeing can range 
from O'/4 to >1/0 in the NIR). 

In order to overcome these limitations, we present 
a new tool named GalPaK 3D (Galaxy Parameters and 
Kinematics) designed to be able to disentangle the 
galaxy kinematics from resolution effects over a wide 
range of conditions. This is not the first code to model 
galaxy kinematics from 3D data (e.g. the TiRiFiC 
package, which perform s tilted ring model fits to three- 
dimensional radio data; |Jozsa et ah 2007"), but this code 
performs disk model fits to three-dimensional IFU data 
cubes, for the first time,[^| whereas all other modeling of 
IFU-data so far h ave worked from the two-dimen sional 
velocity field (e.g. Cresci et al.||2009| Epinat et al. 


2009; 

Davis 


Davies et al.| 20111 |Andersen & Bershady||2013| Davis 

et al.||2013|). 

'Phis paper is organized as follows: we describe the 
GalPaK 315 algorithm in Section [2| We present some test 
case examples in Section [3] We present results from an 


4 Available a t http: //galpak. irap. omp. eu/ 

5 |Law et al. [ (|2012|) made an attempt at 3D fitting, albeit not 
self-consistently. 
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extensive analysis of 1728 synthetic galaxies in Section |4j 
where we discuss the impact of the accuracy in the point- 
spread-function (PSF) characterization. In Section [5j 
we present an analysis of data cubes generated from hy- 
drodynamical simulations of isolated disks from Michel- 
Dansac et al. (in prep.). We summarize this paper in 
Section[6j Throughout, we use the following cosmological 
parameters: H 0 = 70 km s _1 , = 0.7 and = 0.3. 

2. THE GALPAK 30 ALGORITHM 

In this section we outline the algorithm principles, 
which are designed to be able to determine galaxy 
morphokinematic parameters from the three-dimensional 
data cube directly. We discuss the merits of using the 
parametric forward fit and its limitations. 


2.1. A Parametric Galaxy Model in Three Dimensions 

Traditionally, kinematic analyses use two-dimensional 
maps generated by applying line-fitting codes to deter¬ 
mine the line wavelength centroids and widths, which are 
only considered to be reliable for spaxels with sufficiently 
high signal-to-noise (S/N) ratios. This S/N condition is 
easily met at low redshifts, but is harder to meet for 
small, high-redshift galaxies. In principle, the choice to 
work in 2D or 3D space is equivalent, but we will show 
that our method can work in the regime (on the spax¬ 
els) where the signal-to-noise ratio (SNR) per pixel (SNR 
pixel” 1 ) is not sufficient for line-fitting codes, which re¬ 
quire a minimum SNR on all spaxels. 

When the PSF FWHM can be characterized to suffi¬ 
cient accuracy^] (within 10% or 20%; see Section]^]), one 
can take its characteristics, together with the instrumen¬ 
tal line-spread function (LSF), into consideration and re¬ 
cover the intrinsic modeled galaxy parameters. The al¬ 
gorithm uses the spectral LSF and the spatial PSF to 
generate a three-dimensional kernel whose characteris¬ 
tics are set for the given instrument (or a user-generated 
instrument module). 

While a full deconvolution of hyperspectral cubes 
would be preferred, it is usually a challenge mathemati¬ 


cally (a new meth od has b een proposed recently by Vih 
leneuve fe Carfantan|2014), and a forward convolution of 


a parametric model offers a very useful alternative. This 
forward convolution gives us the opportunity to estimate 
intrinsic modeled kinematic parameters in a wider range 


Bouche et al.|| 2013; Peroux et al. 2014| 

Schroetter et al. 

2015 Martin & Soto 2015 for first app 

ications). 


model, and we focus here on a galaxy disk model for 
emission-line surveys, but the algorithm is adaptable to 
other situations. In order to construct a modeled galaxy 
in the observational coordinate systems ( x , y, A), we start 
by generating a three-dimensional galaxy model in a Eu¬ 
clidian coordinate system (x, y , z), where the z-axis is 
normal to the galaxy plane ( x , y). We apply a radial flux 
profile /(r), from one of the traditional Gaussian, expo¬ 
nential, and de Vaucouleur choices as parameterized by 


the Sersic (1963) profile: 


I{r) = I e exp (r/R e ) 1/n - 1 ) 


(1) 


6 The PSF shape matters mor e th an the level of accuracy on the 
FWHM, as discussed in Section 14.41 


with n = 0.5, 1.0, and 4.0, respectively, where R e is the 
effictive radius, /tot the total flux, and b n such that R e 
is equivalent to the half-light radius i?i/ 2 , and I e the SB 
at R e . For n = 0.5, 1.0, and 4.0, the constant b n is 
0.69, 1.68, and 7.67, respectively, from b n ~ 1.9992 n — 
0.3271. The Sersic index n is kept fixed given the large 
degeneracies it creates with other parameters, such as 
the galaxy half-light radius. This degenaracy is due to 
the fact that the SB profiles around R e are c lose to on e 
another fo n = 0.5, 1.0, or 4.0 as noted in Graham et al.| 
([20051. 

To this two-dimensional disk model, we add a disk 
thickness h z . We adopt a Gaussian luminosity distribu¬ 
tion perpendicular to the plane, I(z) tx exp(—z 2 /2 h z ), 
defining h z as the characteristic thickness of the disk. 
GalPaK 3D also allows the user to choose an exponen¬ 
tial I(z ) a exp(— \z\/h z ) or a sech 2 distribution I(z) oc 
sech 2 {z/h z ). We set the disk thickness to h z = 0.15i?i/ 2 
where i?i/ 2 is the disk half-light radius. This choice cor¬ 
responds to h z ~ 1 k pc, ty pical of high-redshift edge- 
on/chain galaxies (Elmegreen & Elmegreen|j2006). At 
this stage, we have a disk model in Euclidean coordi¬ 
nates that accounts for the flux distribution only. 

For the gas kinematics, we create three kinematic cubes 
in the same spatial coordinate reference frame for the 
velocities v = (V x ,V y ,V z ) assuming circular orbits. The 
rotational velocity vyr) with a maximum rotation veloc¬ 
ity I'm ax can have several fu nction al for ms: it can be 
an arctan velocity profile (e.g . Puech et al. 2008), an in¬ 
verted exponenti al (|Feng fe Gallo 20 lTb, or a hyperbolic 
tanh profile (e.g. Andersen &; Bershady 20131 : 


v(r) = V max — arctan (r/r t ) 
7r 

v{r) = V m 
v{r) = V m 


[1 - exp(r/r t )] 
tanh(r/r t ) ‘tanh 


‘arctan’ 

‘exp’ 


( 2 ) 

( 3 ) 

( 4 ) 


where r is the radius in the galaxy x, y plane, r t is the 
turnover radius, and IA max is the maximum circular ve¬ 
locity. These choic es are more extensively discussed in 
Epinat et al. (2010), but it is worth noting that the ‘exp’ 
and hyperbolic rofation curves have a sharper transition 
around the turnover radius. We stress that our parame¬ 
ter Vmax is not the projected asymptotic velocity, but is 
the true asymptotic velocity irrespective of the inclina¬ 
tion. 

Another option, called “mass,” assumes a constant 
light-to-mass ratio and sets v(r) from the enclosed 
light/mass /(< r) profile 


v{r) oc y /(< ^ ‘mass’ (5) 

where r is the radius in the galaxy x,y plane and V max 
normalizes the profile. This option has a rotation curve 
that peaks at some radius (set by the half-light radius), 
decreases at larger radii, and is to be preferred for nu¬ 
clear disks or when there is no significant dark matter 
component. 

We then rotate the disk model around two axes ac¬ 
cording to an inclination (*) and position angle (PA, an¬ 
ticlockwise from y) and create a cube in x, y , and A using 
the three intermediate 2D maps: the flux map, the ve¬ 
locity field, and the dispersion map (otot)- The flux map 
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TABLE 1 

Default Range on Each Parameter 


Parameter 

Min 

Max 

x c 

1/3 Npix x 

2/3 TVpixa; 

Vc 

1/3 Apixy 

2/3 A/pixy 

Zc 

1/3 Afpix z 

2/3 Npix z 

Flux 

0 

3 X ^ij,k( v i,j,k) 

Pl/2 

0.2 spaxel 

4 " 

Inch (deg) 

0 

90 

PA (deg) 

-180 

180 

n 

0.01 spaxel 

1 " 

Vmax (km S -1 ) 

-350 

350 

( T 0 ( kms -1 ) 

0 

180 


is obtained from the rotated flux cube summed along the 
wavelength axis. The velocity field is obtained from the 
flux-weighted mean V z velocity cube. The total (line-of- 
sight) velocity dispersion cr tot is obtained from the sum 
of three terms (added in quadrature). It includes (i) the 
local isotropic velocity dispersion <Td driven by the disk 
self-gravity, which is oyi (f)//i -x = V(r)/r for a compact 
thick” or large “thin” disk flGenzel et al.||2008 Binney 
i Tremaine! 20081 [Davies et ah 2011)V (ii) a mixing term 


& 


a m , arising from mixing the velocities along the line of 
sight for a geometrically thick disk, which is obtained 
from the flux-weighted variance of the cube V Zl and (iii) 
an intrinsic dispersion ( a 0 ) —which is assumed to be 
isotropic and constant spatially— to account for the fact 
that high-redshift disks are dynamically hotter than the 
self-gravity expectation. Indeed, this turbulence term cr 0 
is often observed to be ~ SO-SOkms -1 in z > 1 disks 
(Law et al. 12007 120121 iGenzel et al. 1120081 iCresci et al.l 
20091 |F5rster Scnreiber et~ al.||2UU91 1 Wright et al.[ 2009; 
Epinat et al.|2010[ |2012| [Wisnioski et ai.||201ip and thus 
dominates the other two terms since the mixing term a m 
is typically ~15 kms -1 and the self-gravity term is 
typically 10-30kms _1 . 

To summarize, the flux profile can be chosen to be ‘ex¬ 
ponential’ (n = 1.0),‘gaussian’ (n = 0.5), and ‘de Vau- 
couleur’ (n = 4.0); the velocity profile v(r) can be arctan 
(“arctan”), inverted exponential (“exponential”), hyper¬ 
bolic (“tanh”) or that of mass profile (“mass”); and the 
local dispersion can be that of the thin or thick disk. 
There are in total 10 free parameters^] to be determined 
from the data. The 10 parameters are the x c , y c , z c posi¬ 
tions, the disk half-light radius i?i/ 2 , the total flux f to t, 
the inclination i, position angle PA, the turnover radius 
rt, the maximum circular velocity V max , and the one¬ 
dimensional intrinsic dispersion o 0 . We will refer to the 
last two (Vmax, a 0 ) as kinematic parameters. Finally, the 
simulated galaxy is convolved (in 3D) with the PSF and 
the instrumental LSF specific for each instrument]^] The 
3D convolution is performed using fast Fourier transform 
(FFT) libraries. 

2.2. The Markov Chain Monte Carlo (MCMC) 
Algorithm 

In order to determine the 10 free parameters on hy- 
perspectral cubes, one needs an algorithm that is inde¬ 
pendent of initial guesses on the parameters, that can 

7 There are only nine free parameters when the “mass” profile 
is used for v[r) since the turnover radius rt is irrelevant. 

8 The user can choose a Gaussian PSF, a Moffat PSF. The PSF 
can be circular or elliptical with a user-defined axis ratio b/a. 


converge even in the presence of local minima, and that 
can handle low S/N data. This is particularly difficult 
for traditional minimization methods because the y 2 hy¬ 
persurface is very flat (outside the shallow well near the 
optimum parameters), and as a result the minimization 
algorithm tends to not converge and be very susceptible 
to local minima. 

Here we use an algorithm to optimize the parameters 
using Bayesian statistics with flat priors on bound inter¬ 
vals for each of the parameters. The algorithm constr ucts 


MCMCs with a Metropolis-Hasting (MH) sampler (Hast¬ 
ings fl 9 70). At each iteration we compute the new set of 


parameters Xi+\ from the last Xi set with a proposal dis¬ 
tribution P from which to draw: 


x i+ i = sti + hP(x i+ i\xi), 


( 6 ) 


where the new set of parameters is accepted or rejected 
as in any MH algorithm. The new proposal set of param¬ 
eters Xi + 1 is then accepted or rejected according to the 
posterior distribution, which amounts to the likelihood 
£ oc exp — x 2 in the considered case of flat priors on the 
parameters. In other words, we assume that the pixels 
are independent and that noise properties are Gaussian, 
which is appropriate for optical/NIR data taken in the 
background-limited regime, and the user can provide the 
full variance cube. More appropriate likelihood functions 
for low counts with Poisson noise can be found in[Migliell| 
(1999). 


The scaling vector h in Equation kfi is derived from the 
variance on the flat (uniform) prior"distributions, whose 
boundaries are adjustable (the default values are listed 
in Table [lj. The user may need to rescale the vector h 
in order to have acceptance rates between 20% and 50%. 
Convergence is usually achieved in a few hundred to a 
few thousand iterations, even though we typically let the 
algorithm run for 15,000 iterations. 

In principle, one has the fr eedom to use any proposal 
distribution P (e.g. MacKay 2003). A Markov chain is 
said to converge to a single invariant distribution (the 
posterior probability) when the state of the chain per- 
sits once it is reached and is said to be ergodic when the 
probabilities x n converge to that invariant distrib ution 
as n —> oo, irrespectively of the initial parameters (Neal 
1993) ®] In addition, if the sampler satisfies the following 
conditions P(x\x') = P{x'\x ), as we have used, the algo¬ 
rithm reverts to the Metropolis method, which satisfies 
the two conditions. In practice, however, one also needs 
a distribution that probes the parameter space efficiently 
in order to converge in a reasonable number of cpu hours, 
regardless of the initial parameters. 

A common proposal distribution is the uniform dis¬ 
tribution that gives equal probabilities to all possible 
values. The Gaussian proposal distribution P(x'\x) = 
Af(x , 1) is probably the most commonly used and is pop¬ 
ular but has one major drawback: the Gaussian distri¬ 
bution is rather narrow such that the algorithm becomes 
sensitive to the initial conditions, making the time to 
convergence to the optimum values very sensitive to the 
initial guess. If the width of the proposal distribution is 
small, the convergence is too slow/large, and when it is 


9 Available 
res-mcmc.html 


at 


http://www.cs.toronto.edu/~radford/ 
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large (for convergence purposes), it will lead to low ac¬ 
ceptance rate and poor efficiencies for convergence. To 
remedy this problem, one could use a mixed distribution 
with a Gaussian draw, say, 90% of the time and a uniform 
draw 10 % of the time, allowing the chain to escape from 
a local minimum. Compared to the Gaussian proposal, 
the mixed distribution has one additional parameter that 
needs to be fine-tuned to the problem, such as the mixing 
ratio. _ 

( 119871 ), 


A third option, as advocated by Szu & Hartley 


is to use a draw from a Cauchy distribution that has 
by definition longer wings (i.e. P is a Lorentzian profile 
where P(x'\x) oc 7 2 /[ 7 2 + ( x' — a;) 2 ])- The Lorentzian 
wings are important, allowing the chain to make large 
jumps during the initial “burn-in” phase and ensuring 
rapid convergence of the chain with no sensitivity to the 
initial parameters. Another advantage of a Cauchy pro¬ 
posal distribution is that it has only one parameter, 7 , 
compared to the mixed one. 

We tested these various choices on simulated cubes and 
found that the Cauchy proposal distribution converged 
faster than the other methods and was least sensitive to 
the initial parameters. In other words, with the Cauchy 
nontraditional proposal distribution, a few hundred to a 
few thousand steps of the MCMC are required to pass 
the burn-in phase depending on the S/N of the data, and 
it is the user’s responsibility to confirm that the MCMC 
chain has converged. Thus, we typically run the chain 
through 10,000 or 15,000 steps to robustly sample the 
posterior probability distribution. 

The “best-fit” values of the parameters are determined 
from the posterior distributions. We use the median and 
the standard deviation of the last fraction (default 60%) 
of the MCMC chain to determine the ‘best-fit’ param¬ 
eters and their errors, respectively. One can also use a 
fraction (default 60%) of the MCMC chain around the 
minimum x 2 - The full MCMC chain is saved such that 
the user can use his/her preferred technique. 

The algorithm is implemented in Python and uses the 
standard numpv an d sci py librai res. In addition, it uses 
the bottl eneck P (Frig o fe Johnson 2005) and FFTw p] 
libraries (|Frigo & Johnson||20I2f iiToraer to speed up 
certain matrix operations and the PSF+LSF convolu¬ 
tion, respectively. It requires FITS files as inputs. The 
algorithm is modular so that the user can add specifica¬ 
tions for other instruments. The online documentation 
describes the syntax, and it takes about 2, 5, and 10 min¬ 
utes on a laptop (at 2.1 GHz) to run 10,000 iterations on 
a data cube with 30 3 pixels, 40 3 pixels, and 60 3 pixels, 
respectively. In other words, the computation time scales 
as t oc lVpi x log(iVpi x ) where N plx is the number of pixels, 
showing that the FFT calculation dominates. 


3. HIGHLIGHT APPLICATIONS 
3.1. Example on 2D data 

Before applying the tool on 3D data, it is important 
to validate the method on simpler data sets, such as 
two-dimensional imaging data. We thus wrote a two- 
dimensional version of the algorithm, GalFit 2D , one that 
does not include the kinematic, which is in essence sim- 


10 Available at https://pypi .python. org/pypi/Bottleneck 

11 Available at https://pypi.python.org/pypi/pyFFTW 


ilar to other par ametric algorithms (e.g. Simard 1998 
Peng et ap|2002 ), apart from the Bayesian approacKT 


"Figure [H shows a comparison between the derived mor¬ 
phological parameters from two data sets of very dif¬ 
ferent resolution. Panel (a) shows a Canada France 
Hawaii Telescope (CFHT) /-band im age of the z ~ 
0.2 g alaxy SDSSJ165931.92+023021.92 ( |Kacprzak et~ah 


2014|. Panel (e) shows an r-band image of the galaxy 


from the Sloan Digital Sky Survey (SDSS) at a spatial 
resolution of T/l. For each data set, we show the fit¬ 
ted (convolved) model, the residual map, and the one- 
dimentional SB profile. One sees that the intrinsic mod¬ 
eled morphological parameters found from the SDSS 
data (PSF FWHM=l7l) are in good agreement with the 
higher-resolution data (PSF FWHM=0'/7). Moreover, 
the residuals in both data sets show the spiral arms and 
a minor merger (or a large clump) in the southern part of 
the galaxy, showing that a smooth axis-symmetric model 
can be used to unveil asymmetric features. 


3.2. Example on a mock cube 

Figure [2] shows an example of a mock disk model with 
a low SNR (SNR pixel -1 of 4 in the central pixel) drawn 
from the set presented in § [4] and generated at 170 reso¬ 
lution. The top, middle, and bottom rows show the flux 
map, the velocity map, and the apparent velocity profile 
V z (r) across the major axis, respectively. From left to 
right, the panel columns show the data, the convolved 
model, the modeled disk (free from the PSF), and the 
high-S/N high-resolution reference data (PSF=0715 and 
S/N=100). In the bottom panels, the solid red curves 
correspond to the reference rotation curve (obtained from 
the reference data set), and the triangles represent the 
apparent rotation curve. These rotation curves show that 
the recovered kinematics from the modeled disk (intrin¬ 
sic or unconvolved model) shown in the third column is 
in good agreement with the reference data (last column) 
in spite of the low spatial resolution (170) and the low 
SNR in the mock data set. 

This synthetic data cube was generated with a flux 
profile with Sersic index n = 1 and half-light radius 
R 1/2 = 075, corresponding to 2.5 MUSE/KMOS pixels), 
an “arctan” velocity profile with U max = 200 kms -1 , a 
thick disk with a velocity dispersion a 0 = 80 kms - , an 
inclination i = 60°, a PA= 130°, and with instrumental 
specifications for the new VLT MUSE instrument (072 
pixel -1 , 1.25Apixel -1 , LSF=2.14 pixels). The integrated 
total flux is 10 -16 ergs -1 cm -2 , and the synthetic noise 
per pixel is a = 5 x 10 -2 ° ergs -1 cm -2 A -1 . 

The synthetic data cube is also displayed in Fig¬ 
ure [3j which shows three one-dimensional spectra (a) 
taken at the three locations labeled in the image 
shown in panel (b). Panel (c) shows a 3D rep¬ 
resentation of the data (blue) with the model over¬ 
laid (red) made with the “visit” software j 12 ] where the 
light/dark areas corresponds to two cuts at fluxes of 6 
and 8 xlO -20 ergs -1 cm -2 A -1 , i.e. an S/N pixel -1 of 
1.2 and 1 . 8 , respectively. 

We ran the algorithm with 15,000 iterations, and Fig- 
ure[4]shows the MCMC chains for the 10 free parameters 
along with the y 2 evolution in the bottom panel. The 

12 Available at http://visit.llnl.org/ 
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Fig. 1. — Application of the two-dimensional version of the MCMC algorithm (‘GalFit2D’) on the z ~ 0.2 SDSSJ165931.92+023021.9 with 
nn r =18.40 mag from Kacprzak et al. (2014). Similarly to GalPaK 3D , Galfit2D performs a parametric fit with an MCMC algorithm using 
set surface brightness profiles convolved with the seeing. The top row shows the result from archival CFHT /-band taken at a resolution 
of 0'/7. The bottom row shows the result from the SDSS r-band image that has a resolution of l'/l. Panels (a) and (e) show the data. 
Using an exponential profile, panels (b) and (f) show the seeing-convolved model; (c) and (g) the residuals, i.e. data-model normalized to 
the pixel noise rj. and (d) & (h) the one-dimensional SB profile. The recovered intrinsic disk scale length R f \ is about 1” in both cases, in 
spite of the different spatial resolution. 
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Fig. 2.— Example of the algorithm application on a disk model simulated with a seeing of 1//0 (FWHM) and a flux of 10 16 ergs 1 cm 2 , 
and an S/N pixel -1 of ~4 at the brightest pixel. The top, middle, and bottom rows show the flux map, the velocity map, and the apparent 
velocity profile V z (r) across the major axis, respectively. From left to right, the panel columns show the data, the convolved model, the 
modeled disk (free from the PSF), and the high-S/N high-resolution reference data (PSF^O'/lh and SNR=100). In the bottom panels, 
the solid red curves correspond to the reference case, the triangles represent the apparent rotation curve, and the dotted lines show the 
apparent V max sin i. One sees that the velocity profile from the modeled disk (third column) is in good agreement with the reference data. 
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(a) 



(b) 



(c) 

Fig. 3.— For the synthetic data in Figure [2] we show three 1- 
dimensional profiles (panels (a)) comparing the data (thin line) 
and the model (thick line) taken at the location labeled “1,”, “2,” 
and “3” in panel (b). Panel (c) shows a 3D representation of the 
data (blue) with the model overlaid (red) where we used two flux 
levels at 5xlO -20 and 1.5xl0 -19 ergs -1 cm -2 A -1 . The cube 
orientation is shown, where the wavelength axis is the ^-direction. 
(An interactive version of this figure is available in the published 
online version) 


values of the fitted parameters (and their errors) shown 
by the black lines (gray lines) are computed from the 
median (standard deviation) of the last 60% iterations 
of the posterior distributions. The recovered parameters 
are listed in Table [2] and show good agreement between 
the input and recovered values. 


TABLE 2 

Comparison between the model input values and the 
RECOVERED VALUES WITH 1 a ERRORS AND CONFIDENCE INTERVALS 
(Cl) FOR THE EXAMPLE SHOWN IN FIGURE [2] 


Parameter 

Input 

Output 

[95% Cl] 

x c (pixel) 

15 

15.05i0.09 [14.87; 15.24] 

y c (pixel) 

15 

15.06i0.09 [14.89;15.23] 

z c (pixel) 

15 

15.05i0.07 [14.92;15.19] 

Flux (10 -16 ) 

1 

1.06i0.03 

[1.01; 1.09] 

/r, (arcsec) 

0.82 

0.85i0.04 

[0.78;0.95] 

Inch (deg) 

60 

62i3 

[58;68] 

PA. (deg) 

130 

126i2 

[123; 130] 

r t (pixel) 

1.35 

1.32i0.42 

[0.8;2.47] 

Umax (kms -1 ) 

200 

202i22 

[172;257] 

g q (kms -1 ) 

80 

82i5 

[73;90] 


Figure [5] shows the joint distributions for the radius, 
PA, inclination, maximum velocity, and dispersion pa¬ 
rameters. The estimated parameters and their respective 
1 er error are shown as a solid line and dashed line, respec¬ 
tively. This figure shows a clear covariance between the 
turnover radius and the asymptotic velocity Vm ax , and a 
small covariance between the inclination and U max . The 
users of GalPaK 3D are strongly advised to confirm the 
convergence of the parameters using diagnostics similar 
to Figure [4] and to investigate possible covariance in the 
parameters, as these tend to be data specific, using di¬ 
agnostics similar to Figure [5} 

4. TESTS WITH MOCK DATA CUBES 

In order to characterize the performances and limita¬ 
tions of the GalPaK 3 # algorithm statistically, we gen¬ 
erated a set of 1728 cubes again with a MUSE configu¬ 
ration over a grid of parameters listed in Table [5] The 
synthetic cubes were generated with noise typical to a 
1 hr exposure with MUSE corresponding to a pixel noise 
of er = 5 x 10 20 ergs -1 cm -2 A -1 . We use a range of 

TABLE 3 

Range of parameters for the 1728 mock galaxies 


Parameter Grid Values 


Flux (10 17 ergs 1 cm 2 ) 

3, 6, 10, 30 

Seeing (”) 

0.6, 0.8, 1.0, 1.2 

Redshift 

0.6, 0.9, 1.2 

# 1/2 (kpc) 

2.5, 5, and 7.f a 

#1/2 (”) 

073, 076, and 17( n 

Inch (deg) 

20, 40, 60, 80 

PA (deg) 

130 

rt (”) 

o.i-o.:£ 

Umax (kms -1 ) 

110, 200, 280 

G 0 (kms -1 ) 

20, 50, 80 


a Exact va lue to satisfy the size-velocity scaling relation |Dutton 
|et al.|2011 ). 

b Exact value will depend on the redshift. 

c Exact value to satisfy the scaling relation between the galaxy 
size and the inner gradient (lAmorisco Sz Bertin 2010) using rt = 
jRd/1.8. 
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Fig. 4.— Full MCMC chain for 15,000 iterations for the example shown in Figure [2] Each of the small panels corresponds to one 
parameter. One sees that the “burn-in” region is confined to the first 1000 iterations. The estimated parameters are shown with the red 
line and are calculated from the last 60% of the chain. The gray lines show the 1 cr standard deviations, and the dotted lines show the 95% 
confidence interval. Note the fitted flux value is (1.06 =b 0.03) x 10 -16 ergs -1 cm -2 , which is found from the sum of the pixel values (here 
~ 8.5 10 —17 ergs -1 cm -2 A -1 ) times the 1.25 A per spectral pixel. The bottom panel shows the x 2 evolution relative to the minimum, 
l°g[x 2 — XmiiJ' We use this non-standard metric in order to show that the variations of the % 2 around the minimum which are 3 to 4 orders 
of magnitude smaller, reflecting a very flat hypersurface. Hence, a plot of % 2 or of the likelihood would show a straight line. 



velocity disoersion 

Fig. 5.— Joint distributions for the radius, PA, inclination, maximum velocity, and dispersion parameters for the example shown in 
Figure [2][4| The estimated parameters and their respective 1 a error are shown as a solid line and dashed line, respectively. One sees 
that the traditional degeneracy between Vhiax and the inclination i is broken, thanks to our 3D method, but the seeing leaves a significant 
correlation between Vhiax and the turnover radius rt. The presence of this degenerancy is data specific and seeing specific, not a generic 
feature of the algorithm. 
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inclinations i from 20° to 80°. We use a range of disk 
sizes, with half-light radii R 1/2 = 0"3, 076, and 170 cor¬ 


range of observed sizes at 2 ~ 1 (e. 

g. Trujillo et al. 2006 

Williams et al. [2010} Dutton et al. 

2011). 


ing relation (Equation 8 of|Dutt on et al. 2011) and its 
redshift evolution (Equation 5 ot |L)utton et al. 2011) to 
set the rotation kinematics (Vm a xjl In particular, the 
sizes R 1/2 = 2.5, 5, and 7.5 kpc correspond to U max val¬ 
ues ranging from ~ 100 to 250 kms -1 . We use “arctan” 
rotation curves to generate our mock data cubes, and 
we have verified that our results remain the same with 
“exponential” rotation curves. 

We use the scaling relation between the turnover ra¬ 
dius r t and the disk scal e length /fo that exists fo r disk 
galaxies (e.g. Figure 1 of Amorisco & Bertin|20 10) to set 
the turnover radius rt- In particular, we set r t to Ad/1-8 
where the 1.8 factor □ is determined empirically for the 
arctan rotation curve to satisfy the linear correlation be¬ 
tween the galaxy disk scale-length = R 1 / 2 / 1-68 and 
Rq, define d as t he radiu s r where V(r) = 2/3 Vmax 
(Amorisco & Bertin||2010). 


For each of the galaxy sizes, the disk thickness is h z = 
0.15 i?i /2 , he. ranging from 0.4 to 1.3 kpc, bracketing 
the average values of h r ~ 1 kpc, fo und for high-redshift 
edge-on/chain galaxies (Elmegreen & Elmegreen 2006). 

We used fluxes for an [OllJ (A3727) emission line, 
expected to lie in the MUSE spectral range at red- 
shifts between 0.6 and 1.2, with integrated fluxes from 
3 x 10” 17 ergs” 1 cm -2 to 3 x 10” 16 ergs” 1 cm” 2 cor- 


responding to the range of observed v alues (e.g. Ba- 
|con et aT 2015| jComparat et al.||2015[ and references 
therein). We use a constant noise value per pixel of a = 
5x 10” 20 ergs” 1 cm” 2 A” 1 , in order to simulate the noise 
level of a 1 hr exposure, but we stress that the algorithm 
accepts variance/noise cubes to account for pixel-to-pixel 
noise variations. In addition, we generated cubes with 
very high S/N (S/N=100, flux= 3 x 10” 15 ergs” 1 cm” 2 ) 
and with a seeing typical of AO conditions, with a PSF 
FWHM of O'/15. These will serve as reference data sets. 

4.1. Surface Brightness and Signal-to-noise Ratio 
One could imagine that the S/N in the recovered pa¬ 
rameters be a function of the average S/N pixel” 1 , or the 
apparent SB since the observed central SB scales directly 
with the SNR in the central pixel. But clearly the com¬ 
pactness of the object with respect to th e seeing plays 
a large role (as discussed in Driver et al. 2005| Epinat 


et al. 2010). Very compact objects (compared to the 
Beam or the PSF) have high SB by definition (and high 
S/N pixel” 1 ), but the morphology and/or kinematic in¬ 
formation may be lost owing to the beam smearing. On 
the other hand, very extended objects have low surface 
brightness (and low S/N pixel” 1 ), but have many pixels 
in the outer regions (with low S/N), where most of the 
information on the galaxy is located and not affected by 
the beam. 

Before illustrating this point, it is important to define 

13 For ‘exponential’ rotation curves, one should set rt to X 0.9 
in order to satisfy the scaling relation; for ‘tanh’ rotation curves, 
one should set rt to Ft ( \ X 1.25. 


commonly used terms such as the SB of galaxies. From 
any light profile I{r) such as given by Equation [lj there 
are many ways to define galaxy SB, such as J e the SB at 
the effective radius R e , I 0 the intrinsic SB at the central 
pixel, A 0 the observed SB at the central pixel, and SB^ 
the average SB within the intrinsic half-light radius R 1 / 2 ' 


SB 


l/2,conv — 


0-5 F tot 

7F ^l/2,conv 


(7) 


where F to t is the galaxy total flux. A related quantity to 
Equation [7] is the observed SB, defined as : 


SHl/2,obs — 


0.5 F u 


S 


1/2,obs 


( 8 ) 


where F tot is the galaxy total flux and S 1 / 2 ,obs the galaxy 
apparent area given by = ira b where a and b are the 
observed major and minor semiaxes, respectively, of the 
galaxy. The relations between these various definitions 
are described in the appendix. 

To illustrate the point made at the beginning of this 
section, we show in Figurepla) the relative errors Sp/p = 
(Pfit — Pin)/ Pin on some ofthe estimated parameters for 
our mock data cubes generated in Section |d] as a function 
of central SB, SBi/^obs (defined in Equation [8]) . Each 
row shows the relative errors for the maximum circular 
velocity Umax) the size R 1 / 2 , the PA, and inclination i 
from top to bottom, respectively. The crosses, squares 
and circles represent the three subsamples with sizes ~ 
2.5, 5, and 7.5 kpc, respectively. One sees that the errors 
in the morphological parameters (size, PA, inclination) 
do increase toward low SBs, but the threshold point at 
which the relative errors reach ^100% depends on the 
galaxy size, represented by the symbols. This illustrates 
the well-known fact that very extended objects have low 
surface brightness (and low S/N pixel” 1 ) but have many 
pixels in the outer regions that contain useful informa¬ 
tion. 

As argued at the beginning of this section and demon¬ 
strated in Figure |6j SB alone might not be sufficient to 
determine the S/N in the fitted parameters, but the com¬ 
pactness of the galaxy with respect to the beam also plays 
an important role. In Figure pi b), we show the relative 
error Sp/p with respect to the observed SB^^bs times 
the size-to-PSF ratio {R 1 / 2 /Rpsf) 0 ■ The symbols cor¬ 
respond to galaxy subsamples with various sizes as in 
Figure [6](a). The index a was found to be empirically 
~ 1 in order to have the relative errors for each of the 
subsamples follow a similar trend and may differ sightly 
for each of the parameters p. In fact, we find that a 
is approximately 0.8, 1.2, and 1.4 for the size, PA, and 
inclination parameter, respectively. 

These empirical results can be explained by the follow¬ 
ing arguments. The apparent SB within the half-light ra¬ 
dius SBx j 2 conv (Equation [ 7 ]) and the observed SBx^obs 
(Equation |8|) are proportional to the SB (or S/N) of the 
central pixel, A 0 , as shown in the A ppend ix (Equation^. 
In the case of no PSF convolution, Refregier et al. (2012) 
showed that (their Equation 12) the relative error a{a)/a 
on morphological parameters (its major-axis a) scales in¬ 
versely to the central I Q where I 0 is the intrinsic central 
SB (Equation[l][2l). In the presence of a PSF convolution, 
Equation 16 of "Refregier et al. (2012) —which applies 
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here— shows that the relative errors on the major-axis 
a scale as 

— ~ OC Aj 1 (1 + ^PSf/^1/2)' (9) 

where Rpsf is the radius of the PSF (-Rpsf = FWHM/2) 
and i?i /2 the intrinsic half-light radius. 

In our cases, for high-redshift galaxies, the ratio 
I?psf/-Ri /2 is — 1.0 and after performing a Taylor ex¬ 
pansion around -Rpsf/-Ri /2 ~ (1 — %) with x = {R 1/2 — 
-Rpsf)/-Ri /2 and M << V one finds that the factor (1 + 
Rpsf/R\/ 2 ) 4S approximately ~ 2 (1— 2 ;) ~ 2 -Rpsf/-Ri/ 2 - 
Hence, Equation [9] on the errors in the major-axis a be¬ 
comes in the regime where Rpsf/R 1/2 — 1-0: 


er(a) 

a 


oc 


oc 


^V 1 

-Rpsf °) 


-1 


( 10 ) 


which shows that the quality of the estimated morpho¬ 
logical parameters will depend on both the pixel S/N 
(or SB) and the galaxy compactness with respect to the 
beam, -Ri/ 2 -Rpsf, as shown in Figure J6|b) 

In both Figure [6]( a) and [61(b), the gray solid lines show 
the expect ed behavior for the morphological parameters 
(Equation [l0]) and one sees that they agree better with 
the mock aata in the right panels for th e mor ph olog¬ 
ical p arameters. This shows that the Refregi er et al.| 
(2012) formalism describes the relative errors on the mor- 
phological parameters (size, PA, and inclination) rela- 
tivel y w ell, as a first approximation. We note that Equa¬ 
tion [To] is only an approximation to Equation [9] when 
-Ri/ 2 /"Rpsf — 1 and that there might be other dependen¬ 
cies for the other morphological parameters, namely, for 
the PA and for the incli nati on. Her e we refer the reader 
to Table 1 of Refregier et al. (2012) and their Appendix 
for further details; it is beyond th e scope of this paper to 
present a full 3D derivation of the Refregier et al. ([2012 1 
formalism. 

Contrary to the morphological parameters, the errors 
in the kinematic parameter Vmax show strong positive 
(negative) biases in the smallest (largest) mock galaxies, 
represented by the crosses (circles) respectively in the top 
panel of Figure [6j(b) . The positive bias for for the most 
compact galaxies (crosses) with respect to the beam can 
be understood because the Vmax information is located 
mostly in the outer parts of the galaxy, where the S/N 
is too low. The negative bias for the largest galaxies (1” 
in R,\ / 2 ) at low SB is likely due to the spatial cut of our 
mock cubes being too small. 

We will return to the reliability of Vmax in section |4.3| 
and now turn to a more detailed discussion on the reli- 
ability of the parameters (size, inclination, disk velocity 
dispersion, and Vmax)- While we used an arctan rotation 
curve, we note that the following results were found to be 
identical when we used an ‘exponential’ rotation curve. 


4.2. Reliability of morphological parameters 

We have shown in the previous section with Figure [6] 
that the relative errors on the half-lig ht radius f ollow 
appoximately the expectation from the |Refregier et al.| 


(2012) formalism. Here we investigate whether the rela¬ 
tive errors depend on some of the other parameters, such 
as inclination, seeing, and size. 

Figure [7] shows the relative errors (pat — Pin)/Pin for 
several key parameters p. The bottom (top) row shows 
the result for the size parameters R 1/2 (inclination i), 
respectively, as a function of seeing, redshift, inclination, 
and size-to-psf ratio R1/2/ Rpsf- The black curves with 
increasing thickness correspond to subsamples with dif¬ 
ferent SB levels (labeled) where the zero point (dotted 
line) has been offset for clarity purposes. The data points 
represent the median, and the size of the error bars rep¬ 
resent the standard deviation for each of the subsamples, 
where we have typically ~ 100 mock cubes per bin. We 
note that the median standard deviations on the param¬ 
eters (from the posterior distributions) tend to be within 
20 % of these binned standard deviations. 

From this figure, one sees that the GalPaK 3D algorithm 
recovers the intrinsic half-light radius -R 1/2 irrespectively 
of seeing, redshift, and/or intrinsic size. Note that the 
relative errors with respect to size-to-seeing ratio at a 
fixed SB follow roughly the expectation from Equation [9] 
where the factor 1 + (Rpsf/-Ri/ 2) 2 saturates to unity in 
our regime with -R 1 / 2 /-RPSF ~ 1 to 2.5. These results are 
not affected by the choice of the SB profile (Sersic n)[^[ 

From the top row in Figure [7[ one sees that the input 
inclination is recovered except at the two smallest fluxes 
and for the more face-on cases. The reason that the 
algorithm can recover the inclination well is that the al¬ 
gorithm breaks the traditional degeneracy between Vmax 
and i using the SB profile (i.e. the axis ratio b/a ) whereas 
traditional methods fitting the kinematics on velocity 
fields have a strong degeneracy between Vmax and the 
inclination i. 


4.3. Reliability of kinematic parameters 

Figure [8] shows the relative errors Sp/p = (pat— Pin)/Pin 
for the parameters Vm ax (top row) and disk dispersion 
er 0 (bottom row) as a function of seeing, <r 0 , inclina¬ 
tion, and size-to-PSF ratio R1/2/ Rpsf- The curves as 
a function of redshift are not shown, because the relative 
errors do not depend on this parameter as in Figure [7] 
The black curves with increasing thickness correspond to 
subsamples with different SB levels (labeled) where the 
zero point (dotted line) has been offset for clarity pur¬ 
poses. The data points represent the median, and the 
size of the error bars represents the standard deviation 
for each of the subsamples, as in Figure |7[ 

Figure^ top) shows that the GalPaK 3 ® algorithm re¬ 
covers the maximum velocity Vmax irrespectively of see¬ 
ing, disk dispersion, and redshift (now shown) provided 
that the galaxy is not too compact. For small galaxies 
with Ri/ 2 /RpsF less than 1.5, the figure shows that it 
is increasingly difficult to estimate the correct values for 
the most compact galaxies, with large uncertainties and 
significant overestimations of thi s parameter. Th is re¬ 
sult was already pointed out in Epinat et ah j|2010 their 


Figure 13) using 2D kinematic models. Epinat et al. 


(2010) also noted that using a simple flat rotation curve 
to model the disk, the maximum velocity Vmax can be 


14 A curve-of-growth analysis on the two-dimensional flux map 
can sometimes yield a constraint on the Sersic index n and a more 
accurate determination of the intrinsic half-light radius (-# 1 / 2 )- 
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Fig. 6.— Relative errors on the estimated parameters Sp/p, defined as (pnt. — Pin)/Pin- Each row shows dp/p for the maximum circular 
velocity V rr ,ax , the size R1/2, the PA, and inclindation i from top to bottom, respectively. The crosses, squares, and circles represent the 
three subsamples with sizes ~ 2.5, 5, and 7.5 kpc, respectively (Table ^]l. The relative errors in for the morphological parameters (size, 
PA, inclination) are binned. Left/ a): Relative errors as a function of central surface brightness SB 1 /p.obs ' n erg/s/cm 2 /arcsec 2 . Right/ b): 
Relative error as a function of central SB, SB b / 2 , 0 bsi times ( R] / 2 /-Rpsf)°S where R\ j 2 is the galaxy intrinsic half-light radius and Rpsp 
the PSF half-light radius. We found, empirically (see text), that a is approximatel y 0 .8. 1.2, and 1.4 for the size, PA, and inclination 
parameter, respecti vely. These values are close to the expectation of —1.0 of Equation |10| (gray lines) derived for morphological parameters 
in imaging data by |Refregier et al.| (2012). The relative error in U Tn;lx does not follow the expected relation and is subject to strong 
systematics for the smallest and largest mock galaxies (crosses). This is due to V max being constrained in the outer parts of the galaxy, 
where the S/N is thus not sufficient for the compact galaxies or where the mock cube is too small for the largest galaxies. 


recovered with an accuracy better than 25%, even when 
i?i/ 2 /i?psF is less than about ~ 2. 

Figure plbottom) shows the GalPaK 3D algorithm re¬ 
covers theaisk dispersion irrespectively of seeing and red- 
shift (not shown). Given the instrumental resolution of 
MUSE used here (R ~ 130 km/s), small dispersions are 
more difficult to recover. We note that the local disper¬ 
sion is rather sensitive to the instrument LSF FWHM, as 
one might expect. The user can specify more than one 
type of LSF (Gaussian or Moffat), and a user-provided 
vector can be specified if the parametric LSF is not suf¬ 
ficient to describe the instrument LSF. 


4.4. A note regarding trie rbL accuracy 

One could argue that our results are driven by the fact 
that we use the exact same PSF (in 3D) as the one used 
to generate these modeled galaxies. To test the relia¬ 
bility of the algorithm in more realistic situations, when 
the PSF FWHM is not known accurately, we ran the al¬ 
gorithm on the same set of data cubes with a random 
component added to the FWHM of the PSF given by a 
normal distribution with er = 0.1, corresponding to un¬ 
certainties in the FWHM of ~ 20%. We found that the 
accuracy of the spatial kernel (PSF) has little impact on 
the recovered parameters. On the other hand, we find 
that the shape of the PSF is more critical especially for 
the morphological parameter such as the axis ratio b/a 
(or the inclination). We note that sophisticated tools ex¬ 
ist to determine the PSF from faint stars in data cubes 


such as the algorithm of Villeneuve et al. (2011). 


To conclude this section, our algorithm is able to re¬ 
cover the morphological and kinematic parameters from 
synthetic data cubes over a wide range of seeing condi¬ 
tions provided that the galaxy is not too compact and 
has a sufficiently high SB. Thus, for galaxies to be ob¬ 
served with MUSE in the wide-field mode in 1 hr ex¬ 
posure and no AO, we find that the algorithm should 
perform well provided that the SB is greater than a 
few xlCU 1 ' ergs -1 cm -2 arcsec -2 and as long as the 
the size-to-seeing ratio i?i/ 2 /i?psF is larger than 1.5 (or 
A-i/ 2 /FWHM > 0.75). 

5. APPLICATION ON HYDRODYNAMICAL SIMULATIONS 

In the previous section we validated the algorithm on 
synthetic or mock data, which have by definition no de¬ 
fects, i.e. are perfectly regular and symmetric. In order 
to validate the algorithm on more realistic data, we now 
analyze the performance of the algorithm on data cubes 
created from simulated galaxies generated from a hydro- 
dynamical simulation (Michel-Dansac et al., in prep.). 
This is intended to validate the algorithm in the pres¬ 
ence of systematic deviations from the disk model. 

5.1. From Hydrodynamical Simulations to Data cCubes 

The simulation used in this work comes from a set of 
cosmological zoom simulations, each targeting the evo¬ 
lution until redshift 1 of a single halo and its large-scale 
environment. The full sample of simulations is presented 
in details in Michel-Dansac et al., in prep. Here we focus 
on one output of one simulation to complement the test 
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Fig. 7. — Relative error Sp/p —defined as (pg t — Pin)/Pin — for the parameters R \/2 (top) and inclination i (bottom) as a function of 
seeing, redshift, inclination i, and size-to-psf ratio b‘i , 2 /Rpsf (from left to right). The curves with increasing thickness correspond to 
subsamples with different SB levels from SB < 2 X 10 -17 , 2 X 10 -17 < SB < 5 X 10 -17 , 5 X 10 -17 < SB < 1 X 10 —16 , SB > 1 X 10 -16 
ergs -1 cm -2 arcsec -2 respectively, where the zero point (dotted line) has been offset for clarity purposes. SB is the surface brightness 
within the observed half-light radius SB 1 / 2 ,obs times the seeing-to-size ratio 7?i/2/^tpSF> as hi Figure|6j The data points represent the 
median, and the size of the error bars represent the standard deviation for each of the subsample. One sees that the GalPaK 3D algorithm 
recovers the morphological parameters irrespectively of seeing, redshift, and/or intrinsic size. 


cases from section [4] with a more realistic, intermediate 
redshift, star-forming disk galaxy. 

The simulations have been run with the Adaptative 
Mesh Refinement code RAMSES ( Teyssier]|2002 ) using the 
standard zoom-in resimulation technique to model a disk 
galaxy in a cosmological context. Each simulation has 
periodic boundaries and nested levels of refinement in 
a zoom region around the targeted halo, in both DM 
and gas. The refinement strategy is based on the quasi- 
Lagrangian approach. The simulation zooms in a dark 
matter halo inside a 20 h -1 Mpc comoving box, achieving 
a maximum resolution of ~ 200 pc. The virial mass of 
the dark matter halo is approximately 3 x 1O 11 M 0 at 
z = 1, sampled with roughly 600, 000 particles. 

The simulation implements standard prescriptions for 
various physical processes crucial for galaxy formation: 
star formation, metal enrich men t, and k inetic feedba ck 
due to Type II supernovae (Dubois & Teyssier 2008); 
metal advection, metallicity- and density-dependent 
cooling; and U V heating d ue to cosmological ionizing 
background (see|Few et ak[2012j for more details on sim¬ 
ilar simulations but focusing on z = 0 Milky-Way-type 
galaxies). 

The simulated galaxy is a typical z = 1 star-forming 


galaxy with M* = 3 x 10 1 o Mq and a gas fraction of 0.33. 
The galaxy exhibits a disk morphology with spiral arms 
as seen in Figure [9] (top right panel). 

From the outpuTof the hydro-simulation, we generated 
a data cube with the Spectrograph for INtegral Field Ob¬ 
servations in the Near Infrared (SINFONI) instrumental 
resolution and pixel size (0"125 pixel -1 and 2 A pixel -1 ) 
using the star formation rate (SFR) and metallicity in¬ 
formation in each cell. To construct the mock data cube, 
the simulated galaxy is artificially placed at z = 1.3 
(A c ~ 1.5pm for Ha) and rotated with an inclination 
of 60°. Star-forming cells are selected by computing the 
mass of young stars inside each cell of the galaxy. Then, 
we convert this star formation rate into Ha flux using the 


Kennicutt (1998) calibration. For each cell, we also com¬ 
pute the nuxin the [N 11 ] line from the values of the Ha 
flux and the oxygen abund ance following the calibration 
given by Perez-Montero et al. (2009). For each spatial 
element or spaxel, we sum the contribution (to the spec¬ 
trum) of each cell along the line-of-sight. Each contri¬ 
bution has its own line of sight velocity, which blueshifts 
or redshifts the lines. The line width in the spectrum is 
then due to the sum or integral over the cells, which is 
then convolved with the instrumental profile. 
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Fig. 8. — Relative error Sp/p —defined as (pg t — Pin)/Pin — for the kinematic parameters V max (top) and a 0 (bottom) as a function 
of seeing, disk dispersion a 0 , inclination i, and size-to-PSF R| / 2 /Rpsf (from left to right). The curves as a function of redshift are not 
shown because the relative errors do not depend on this parameter as in Figure[7] For the V rn ;\x parameter, the black (red) curves show the 
results when R] / 2 /ffi'SF is less (greater) than 1.5, respectively. One sees that the GalPaK 3D algorithm recovers the kinematic parameters 
irrespectively of seeing and redshift, provided that the galaxy is not too compact with R 1 / 2 /Rpsf larger than 1.5. 


We generated seeing-convolved cubes with seeing of 
0"50, 0765, 0"80, 170 and 172 (corresponding to typi¬ 
cal values in the NIR with SINFONI) and 0715 (cor¬ 
responding to adaptive optic assisted observations) and 
added noise corresponding to a given max S/N pixel -1 . 
Cubes generated with a SNR equal to 100 and a seeing 
of 0715 are used as reference cubes. The final cube size is 
28 x 28 x 30 (in x, y, A directions), but we also produce 
another set of cubes of size 28 x 28 x 200 pixels to allow 
sufficient wavelength baseline for our custom line-fitting 
algorithm that was used to produce the 2D velocity maps 
shown in Figure [9] 

5.2. Application of the algorithm 

Figure |9] shows the results of the GalPaK 313 algorithm 
for a seeing of Of!8 and a minimum SNR pixel -1 of 3 in 
the brightest pixel. As in Figure [2j the top, middle, and 
bottom rows show the flux map, the velocity map and 
the apparent velocity profile V z (r) across the major axis, 
respectively. From left to right, the panel columns show 
the data, the convolved model, the modeled disk (free 
from the PSF), and the high-SNR high-resolution refer¬ 
ence data (PSF=0715 and SNR=100). In the bottom 
panels, the solid red curves correspond to the reference 


rotation curve (obtained from the reference data set), 
and the triangles represent the apparent rotation curve. 
By comparing the two, one sees that the algorithm is 
able to recover the kinematics (third column) in a regime 
where traditional 2D methods (left most column) tend to 
be noisier. In other words, the recovered kinematics from 
the modeled disk (intrinsic or unconvolved model) shown 
in the third column is in good agreement with the ref¬ 
erence data (last column) in spite of the lower spatial 
resolution (078) and the lower S/N in the data set. 

We ran the GalPaK 3D algorithm on the data cubes, 
setting the rotation curve v(r) to an“ arctan” profile and 
setting the Sersic index n to 1.0[//] From the cube with 
a S/N of 100, the inclination found by the GalPaK 3D 
algorithm is 58° ± 2°, and the half-light radius R 1/2 is 
~ 3.4±0.1 kpc (or ~ 074), and its asymptotic maximum 
velocity V max is ~ 215 ± 10 kins -1 , placing it c lose to 
the * ~ 1.5 size-velocity relation of Dutton et al. (2011). 
The asymptotic maximum velocity is close to the one ex¬ 
tracted directly from the simulation, which is 235 kms - . 


15 We also ran the algorithm with “gaussian” profiles with n = 
0.5 leading to very similar results. 
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Fig. 9. — Application of the MCMC algorithm on a disk galaxy generated with the AMR code RAMSES (|Teyssier|[2002| and ‘observed’ 
with a seeing of 0'/8 (FWHM). The maximum SNR is ~3 in the brightest pixel. The top, middle, and bottom rows snow the flux map, 
the velocity map and the apparent velocity profile V z (r) across the major axis, respectively. From left to right, the panel columns show 
the data, the convolved model, the modeled disk (free from the PSF), and the high-S/N high-resolution reference data (PSF=0”15 and 
S/N=100). In the bottom panels, the solid red curves correspond to the reference case, the triangles represent the apparent rotation curve, 
and the dotted lines show the apparent maximum line-of-sight velocity V max sin i. One sees that the the velocity profile from the modeled 
disk (third column) is in good agreement with the reference data (solid line) at ()" 1 5 resolution. 


We repeated the exercise on this simulated galaxy vary¬ 
ing the luminosity (SFR in our case) where the noise 
level is set for a given exposure time corresponding to 
a 2 hr integration with the SINFONI instrument. Fig¬ 
ure 10 shows the maximum signal to noise per pixel (solid 
lines) as a function of the seeing FHWM for five fixed 
SFRs, 5, 10, 15, 30, and 60 M 0 yr _1 , respectively. The 
green region shows the parameter space where the algo¬ 
rithm is able to recover the kinematics parameters within 
20%, from the value determined in the high-S/N cube. 
The yellow region shows the parameter space where the 
algorithm is marginally able to recover the kinematics 
parameters, i.e. within 20%-40% The red region shows 
the parameter space where the algorithm is unable to re¬ 
cover the kinematics parameter, where the relative error 
is larger than 40%. This plot shows that the kinematic 
parameters can be well estimated irrespectively of seeing, 
provided that the SNR is above a critical value (3 in this 
case). Consequently, when the PSF FWHM is slightly 
below the original scientific goal, the optimal observing 
strategy is to integrate longer. 

In the background-limited regime, the S/N per pixel 
scales as oc yft exp , where t e xp is the exposure time. Given 
that the total flux of a circular extended source is F 0 b s ~ 
A 0 7rcr 2 where er 2 = (I? 2 ^ 2 +R.p SF )/1.17 2 , the SNR in the 

central pixel (i.e. the central SB C , or A 0 ) will scale as 

At r 2 

''exp ' pix 


SNR(H 0 ) oc ° 6XP pix = 
Y SBsky t ex p r pi x 


exp 


-Rpsp 


^ 1/2 


' pix 


( 11 ) 


where Rpsf the PSF radius, and R \/2 the object half- 
light and v pjx the pixel size in arcseconds, such that a 
change of (T2 in the PSF FWHM (from O'.'8 to l'/O) cor¬ 
responds to a fraction change of 15% in S/N for a galaxy 
of size R \/2 = O'/6, and accordingly 30% more exposure 
time would be required to reach the same S/N. 


6. CONCLUSIONS 


In this paper we presented an algorithm to constrain 
kinematic parameters of high-redshift disks directly from 
3-dinrensional data cubes. The algorithm uses a para¬ 
metric model and the knowledge of the 3-dinrensional 
kernel to return a 3D modeled galaxy and a data cube 
convolved with the 3D kernel. The parameters are es¬ 
timated using an MCMC approach with nontraditional 
sampling distributions in order to efficiently probe the 
parameter space. 

In summary, 


1. the 2D version of the algorithm is used on an SDSS 
r-band image of a z ~ 0.2 galaxy (Figure [l]) taken 
at l'.'l resolution. We find that the morphology 
is well recovered compared to a higher-resolution 
(O'/7) CFHT image; 

2. using a set of 1728 mock data cubes, Figure [6] shows 
that the accuracy on the recovered parameters de- 
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Seeing (") 

Fig. 10. — S/N pixel -1 as a function of seeing (FWHM) for our 
SINFONI data cube simulated for a 2 hr exposure time. The lines 
correspond to a given SFR at z = 1.3 from SFR=5 to 60 Mq yr —1 . 
The simulated cubes (generated from the hydrodynamical simula¬ 
tion) can be reasonably well fitted by our algorithm provided that 
the S/N pixel -1 (at the central region) is greater than 3, irrespec¬ 
tively of seeing. This diagram applies to galaxies with inclinations 
around ~ 60°. 


pends on the product of the central SB, SB 1 / 2 ,obs 
times the size-to-seeing ratio (Ri/ 2 /Rpsf)~ 1 , fol- 
lowing approximately the analytical expectation of 


Refregier et al. (2012); 


3. from this set of mock data cubes, the morphological 
parameters do not depend on seeing, redshift, or 
the size-to-seeing ratio (Figure [7|; 


4. from this set of mock data cubes, the robustness of 
the algorithm in recovering the kinematics parame¬ 
ters is also independent of seeing and redshift, pro¬ 
vided that the ratio between the galaxy half-light 
radius and the PSF radius (-Ri/ 2 /Rpsf) is larger 
than 1.5 (Figure [8]); 


5. we also find that the accuracy in the recovered pa¬ 
rameters does not depend on the FWHM accuracy, 
but depends more critically on the shape of the 
PSF, except for the disk dispersion a D , which de¬ 
pends critically on the instrument LSF; 


6. using a simulated disk galaxy from the hydro¬ 
simulation of Michel-Dansec et al., which contains 
asymmetric deviations, we found that the kine¬ 
matic parameters can be well estimated irrespec¬ 
tively of seeing, provided that the SNR is above 
a critical value (3 in this case; Figure |To|). Con¬ 
sequently, when the PSF FWHM is sligiitly above 
the original scientific goal (l'.'O instead of Cy/8) the 
optimal strategy is to integrate 30% longer (Equa¬ 


tion 


11) for a galaxy of size i?i/ 2 = 0'.'6. 


In conclusion, the GalPaK 3D algorithm can provide 
reliable constraints on galaxy size, inclination, and kine¬ 
matics over a wide range of seeing and of S/N. However, 
the algorithm should not be used blindly, and we stress 
that users of GalPaK 3D are strongly advised (1) to look 
at the convergence of the parameters (as in Figure [4]); 
(2) to investigate possible covariance in the parameters 
(as in Figure[5]), as these are rather data specific; and (3) 
to adjust the X1CMC algorithm to ensure an acceptance 
rate between 30% and 50%, as discussed in the online 
documentation E3 

Recent a pplications of the Ga l PaK 3D algorithm can 
be found in |Peroux et al.| fl 2013|); IBouche et al.l (|2013|); 
Schroetter et ai.| (|2015|), and|Boiatto ct af. (2015), which 
illustrate the potential in using a global 3D fitting tech¬ 
nique. 
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APPENDIX 

SURFACE BRIGHTNESSES 


For extended sources with total flux Ftot and exponential profiles, i.e. SB(r) = /(r), one can define several measures 
of SBs. We have the following: 


1. The central SB I 0 which is 


In = 


F t( 


2 ^ 


( 1 ) 


in the case of an exponential flux profile /(r) = I 0 exp(— r/Rd) since F to t = / 0 27ri?j, where the half-light radius 
R 1/2 = 1-68 i?d- In the case of a Gaussian flux profile I(r) = I Q exp(—r 2 /2cr 2 ), it is 

Ftot 


In = 


2ircr 2 


( 2 ) 


where the half-light radius R \/2 = 1-17(7. 

2. The average SB within the half-light SB]y 2 : 

cn> _ 0-5 Ftot T 

SB]/2 — -— oc i 0 , 

ttR 2 1/2 

where i?i/ 2 is the true or intrinsic half-light radius (-R 1/2 = 1-68 i?d). 

3. The central pixel SB, SB C : 

SB C = A a 


(3) 


(4) 

where the observed SB profile F 0 \, s (r ) is the convolution of /(r) with the PSF G(r), i.e. F 0 b s (r) ~ 
A 0 exp(—r 2 /2cr 2 ) where now cr contains the contributions from the intrinsic profile and from the PSF via 
(1.17cr) 2 = R\ /2 + R 2 sf (= R 2 1/2 conv ). Rpsf is the radius of the PSF (i? PSF = FWHM/2). 


4. The apparent central SB within the half-light radius -ffi/ 2 , 

0.5 F tot 


SB 


1 / 2 , c 


convi 8 >Hl/ 2 ,conv ■ 

0.5 F tot 


irR 2 

l/2,conv 


r {R\ 


1/2 




(5) 


where -Ri/ 2 , CO nv is the convolved half-light radius. 
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5. the observed central surface brightness within the observed galaxy surface area 5i/2 j0 bsi STb/ 2 . 0 bs : 

_ 0.5 F tot _ 0.5 F tot 

^^1/2,obs 0 7 

*~ ? l/2,obs 

where a and b are the observed major and minor axis, respectively. 


(6) 


The first three (Equations |l}|3| are not observable but can be derived from the total flux F to t and from the galaxy’s 
intrinsic size Rd or R\/ 2 ■ On the other hand, the other two (Equations 4jj5j) are directly observable. 

Naturally, the galaxy apparent area 5j/ 2 . 0 bs is 7Ta b or tt a 2 (5/a); thus, the face-on SB 1 / 2j 
observed SB|/ 2 ,obs (Equation^ are related to one another via the axis ratio fr/qp 7 ] 

From these definitions, we now derive relationships between these variants of SBa 
for intermediate galaxies, the seeing radius i?psF and the galaxy half-light radius R ±/2 are of the same order, i.e. 
-Rpsf/-Ri /2 — 1- Hence, one can write -Rpsf/-Ri /2 — (1 — x) with x = (R 1/2 ~ Rfsf) / R 1/2 and \x\ « 1. 

Since the total flux F tot = 2n a 1 A 0 is also 2n R ^ I 0 , we have 


(Equation |5j) and 
and begin by noting that, typically 


2n lo = A 0 2n 


* 1/2 


R? 


PSF 


In — A n 


2 ln(2) 
1.68 2 i?psF 

R 2 y 


(7) 


which relates the observed S/N in the central pixel A a to the intrinsic central SB I Q . 
The average central surface brightness SB-| / 2iConv within i?i/ 2jC onv is 


SB 


l/2,conv ' 


0.5 F tot 


7r T?i/ 2 1 + (7 ?psf/7?i/2)^ 
0.5 Ftot 1 

l\ 1 2 

-0.25 SB 1/2 x oc I Q x 


R 


1/2 


Fpsf 


7?psf 


( 8 ) 


-0.5 


ln(2) 


which shows that the observed central SB, SB-|/ 2j0 bs: directly maps onto the S/N in the central pixel. 


17 Generally speaking, a = i?i/ 2 ,conv — (Ff/ 2 + ^psf)° 5 an d 
b ~ {R\/2 cos2 (*) + Fp SF ) 0 5 . 













